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ABSTRACT 


The problem of inviscid, steady transonic conical flow formulated 
in terms of the small disturbance theory is studied. The small distur- 
bance equation and similarity rules are presented, and a boundary value 
problem is formulated for the case of a supersonic freestream Mach 
number. The equation for the perturbation potential is solved numer- 
ically using an elliptic finite difference system. The difference equa- 
tions are solved with a point relaxation algorithm that is also capable 
of capturing the shock wave during the iteration procedure by using the 
boundary conditions at the shock. Numerical calculations, for shock 
location, pressure distribution and drag coefficient, are presented 
for a family of nonli.fting conical wings. 

The theory of slender wings is also presented and analytical 
results for pressure and drag coefficients are obtained. 
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CHAPTER 1 


INTRODUCTION 

Transonic flows are characterized by the presence of a region where 
the local flow speed is close to the local speed of sound. In flows of 
this nature shock waves are generally present in the supersonic region 
or at the downstream supersonic-subsonic interface. Mathematically 
this problem is described by a nonlinear equation which is elliptic if 
the flow is subsonic or hyperbolic if the flow is supersonic. 

Intense theoretical and ejcperimental work started in the late 

forties and resulted in a good understanding of the physical phenomenon. 

Newman and Allison^ give a listing of over 300 theoretical papers and 
2 

Cole gives a brief historical survey leading to problems of current 
interest. Those problems are of great engineering interest because 
they occur in a number of practical situations such as flows in nozzles, 
over helicopter rotors, propellers and over airplanes flying close to 
Mach number one. 

The method of the hodograph transformation, which makes the 
equation linear, was the first to yield quantitative information about 
transonic flows over two dimensional bodies. The hodograph method is 
useful in the case of the straight wedge and can also yield a solution 
valid in the vicinity of blunt noses. However, its applicability is 
hindered since boundary conditions, inherently given in the physical 
plane, usually become very complicated in the hodograph plane. The 

3 

limitation is thus in solving direct problems . Garabedian and Korn 
have found this method useful" in computing the inverse problem; the 
airfoil for a given flow field. 


In computing solutions successfully over arbitrary planar shapes 
two basic approaches have been used. Both approaches, the time depen- 
dent techniques and the relaxation methods use finite difference tech- 
niques and are capable of capturing imbedded shocks . The time dependent 
technique (Magnus and Yoshihara)^ integrates the equations of unsteady 
compressible flow forward in time to approach a steady state. This 
method is quite lengthy for general application. A simpler and faster 
method for symmetric bodies was first developed by Murman and Cole^ who 
solved the transonic small disturbance equation for the potential. Sub- 
sequent work (Krupp and Murman, ^ Krupp^) improved this procedure and 
extended it to lifting airfoils. The results have been shown to give 
remarkably good results. 

For three-dimensional transonic flows only a few solutions are 

g 

available thus far. One due to Ballhaus and Bailey, extends the relax- 
ation procedure of Murman and Cole to solve three-dimensional super- 
critical flows over lifting wings with blunt leading edges and non- 
rectangular planforms. However, the solutions available are all for 
subsonic freestream Mach numbers. 

The purpose of this dissertation is to develop on effective 
numerical method for computing three-dimensional transonic flows with 
"supersonic" conditions. More specifically nonlifting conical wings 
with the shock attached to the nose of the wing will be considered. 

This dissertation will be presented in three major parts, consist- 
ing of Chapters 2, 3 and 4. Chapter 2 formulates the boundary value 
problem and gives the basic analytic results needed for the numerical 


computation. Chapter 3 discusses the details of the numerical 
computation and presents the results. Chapter 4 presents the slender 
wing theory and compares the results. 
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CHAPTER 2 


FORMULATION OF GOVERNING EQUATIONS 
2.1 BASIC EQUATIONS 

The basic differential equations are statements concerning the 
conservation of mass, momentum and energy for a fluid, plus an equation 
of state. Since shock waves are present shock jump conditions must be 
added. The shock jump conditions are integral forms of the above con- 
servation laws. Boundary conditions prescribe the flow at upstream 
Infinity and require the flow to be tangent to the body surface. 

The continuity and momentum equations are as follows: 


V.pq = 0 (2.1) 

q-Vq+^Vp=0 (2.2) 


where p is the density, ^ the velocity vector and p the pressure; the 
space coordinates used being non-dimensional. 

The energy equation expresses the fact that total enthalpy is 
conserved along streamlines and it can be shown that it does not jump 
across a shock surface. Assuming the perfect gas law we can write the 
energy equation as. 



(2.3) 


where a = local speed of sound = 



constant 



U = magnitude of the velocity vector at upstream infinity 
a^ = speed of sound at upstream infinity. 

PRECHDiNG PAGE iiLAi^K NOT FlLIviia) 


The only mechanism for entropy production and for introduction 
of rotation in this inviscid model is a shock wave. The entropy, 
though constant along streamlines, does jximp across a shock wave. It 
can be proved that the jump in entropy is of a higher order than the 

9 

equation to be considered. Hence the flow can be assumed isentropic 
and irrotational. 

Using the relation between the speed of sound and the pressure 
we can rewrite the momentum equation as: 

2 


q .Vq = - ^ Vp = - - Vp = - — 

^ ^ p\dp/'^ p 


Vp 


(2.4) 


The continuity equation gives a relation between the density and 
the velocity of the fluid 

V«p^ = ^*Vp + pV.^ = 0 (2.5) 

and multiplying the momentum equation (2.4) by ^ enables p to be 
eliminated using the continuity equation (2.5), giving: 


-^ = a^ 


( 2 . 6 ) 


Introducing the condition for irrotationallty 
V X q = 0 (2.7) 

a velocity potential $(x, y, z) can be defined by the relation 

q = (2.8) 

The components of velocity, in equation (2.6), can be substituted by 
the appropriate partial derivatives of the velocity potential giving a 
scalar equation for the velocity potential. 


(a^ - + (si^ - 

V X/ XX \ y/ yy \ z/ 


)1$ + (a^ - $ + (a^ - $ 

XX \ y/ yy V z/ zz 

-2$x ^y ^xy - 2$x $z ^xz ~ 2$y $yz = 0 


(2.9) 
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If the flow at upstream infinity is uniform with magnitude U and 
parallel to the x axis 

$(x, y, z) = Ux at upstream infinity (2.10) 

On the surface of the body the flow must be tangent to the solid 
surface. This can be expressed as 

q-VS =0 (2.11) 

where S(x, y, z) describes the body surface. 

The equation for the potential, together with the above boundary 
conditions form the basic system describing the flow field. 

2.2 DERIVATION OF THE TRANSONIC SMALL DISTURBANCE EQUATION 

The solution of equation (2.9) gives the exact flow, with the 
assumptions outlined in the preceding section. However, there are no 
practical methods available to solve this equation. We must then look 
for a simpler equation capable of preserving the features of the exact 
equation and which could be solved practically. 

The equation can be simplified by introducing small disturbance 
approximations corresponding to flow over thin bodies of thickness 
26 (6 « 1 ). 

The small disturbance approximation corresponds to an expansion 
procedure, for the velocity potential, valid as certain small param- 
eters approach zero. 

A special procedure is necessary to derive a suitable small- 
disturbance approximation when the flow is transonic as opposed to the 
usual linearized theory. Linearized theory, which has a singularity 
at = 1, cannot be used in the transonic range since the solution it 


gives for the streamwlse velocity perturbation becomes infinite as the 
upstream Mach number approaches one. This violates the assumption of 
small disturbance. 

In the transonic range the expansion should also be as 3 nnptotic 

as the upstream Mach number approaches one, and it will be shown that 

2 

the quantities 6 and M^ - 1 should not go to zero independently. 

The expansion for three dimensional wings of halfspan b and 
thickness 26 (Figure 2.1) has the form: 

$(x,y,z; M^,6) = U |x + £j^(6)^^(x,y, z) + + ... 

y = 3(<S)y (2.12) 

z = 3(6)z 

Linearized theory, where 3(<S) = 1, fails badly near M^ = 1 because 
the orders of certain terms were estimated incorrectly. In order to 
arrive at the correct expansion let us set M^^^ =1. If a suitable 
expansion can be obtained for = 1 then it can be extended in the 
neighborhood of M^ = 1 by considering that M^ ->■ 1 at a certain rate as 
6 -> 0. 

For M^ = 1 this limit process cannot be carried out in the original 
system of coordinates, since this gives the same result as in linear- 
ized theory. Instead rescaled y and z coordinates are used. Linear- 
ized theory Indicates correctly that disturbances spread more rapidly 
in the transverse direction of flow hence we can expect 3(5) 0 as 

- 1 . 

Keeping only the first correction term to the potential and 


dropping the subscript 1 



$ = U{1 + zccl} 

X X 


$ = U e 6 

y y 


$ = U e 3 

z z 


and the energy equation takes the form 


a^ = af - iy-l)}p-zcp + O(e^) 


(2.13) 

(2.14) 

(2.15) 

(2.16) 


Substituting equations (2.13) - (2.16) in equation (2.9), dividing 

2 2 
by U and omitting terms of order e , we obtain 

K-1.- (Y+l)eJ “ (2.17) 


If we now set = 1 in the above equation 


<,y-¥l)z^(4> (p = e3^(0~~ + 

^x XK yy zz 


The distinguished limit occurs if 3 0 as e -»■ 0, in such a way 


that 


2 2 

e3 = e 


or 3 = x/i” 


(2.18) 


The condition that the flow is tangent to the solid surface, 
equation (2.11), fixes the orders of magnitude of 3 and z. 

The wing is described by 


S(x,y,z) = 0 = y - 6f(x, :g) 


(2.19) 

and using equation (2.11) gives the boundary condition for the pertur- 


bation potential on the body. 


e3^(x,0,z) = <Sf^(x, :g-) 


( 2 . 20 ) 
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A comparison of equations (2.18) and (2.20) shows that 


= 6^/3 

fl/3 


( 2 . 21 ) 


6 = 6 ' ( 2 . 22 ) 

for the distinguished limit. This gives the equation for sonic flow 


(Y+1)9? cp = cp-~ + (f) — = 1 

' XX ^yy ^zz • °° 


(2.23) 


The concept of the limit process leading to the expansion must 
now be widened to consider M -> 1 as 6 -»■ 0, this can be expressed as 


= 1 + Kv(6), 


v(6) -> 0 


The quantity K is held fixed in the limit and an order must be found 

for the function v(6) . To do this we go back to equation (2.17) 

2 

rearrange it and substitute Kv(6) for - 1. The distinguished limit 
exists if 

V(6) = 6^''^ (2.24) 


and results in the transonic equation. 

Tk + (y+1)<P 1 (P = 

VI / xj ^xx yy zz 


- 1 


K = 


.2/3 


(the transonic similarity parameter) 


(2.25) 


(2.26) 


The boundary condition on the wing reduces to 


<p~(x.o,5) - tjx, 


(2.27) 


To make the system describing the flow field a complete one we 
must next consider the shock jump conditions. These relations can be 
derived from the surface integral form of equation (2.25) applied 
across the shock surface plus the requirement that the resultant jump 


I 


in velocity Is normal to the shock surface. Using velocity components 
u, V and w to substitute for cp~ and equation (2.25) can be 


written In divergence form. 


u^l + v~ + w~ = 0 
/x ^ ^ 


(2.28) 


The above equation can be Interpreted as the transonic continuity 
equation since it represents the divergence of the mass flux vector. 


With 

Q = Ku - u^^ 1 + vj + wk 

(i,j,k are the unit vectors in the x, y and . z direction respectively) 
equation (2,28) can be written as 

V-Q = 0 

Across the shock surface, Figure 2.2, the jump in the normal 
component of Q is zero and so is the jump in the tangential component 
of velocity. To incorporate these in equation (2.28) let the shock 


surface be given by 


G(x, y, z) = 0 = X - g(y, z) 


Then the unit normal to the shock surface is 


(2.29) 


- VG- - - - 

^ = IWr ~ ^ - 8y J ^ 


[Q • n] = 0 


(2.30) 


j^Ku + u^j + [v]g- + [wlgg = 0 


(2.31) 
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where 


[ ] = jump = ( ) 


( 2 ) 


( ) 


( 1 ) 


Since the tangential component of velocity does not jump across 
the shock 

X n] = 0 (2.32) 

or 


- [v]g~ + [w]g~ = 0 


[ujg^ + [w] = 0 


[u]g~ + [v] = 0 
Substituting equations (2.33) 

[u] + [v] 


V 4 . 2 

Ku + u 


- (2.35) in equation (2.31) gives 

^ + [w]^ = 0 


(2.33) 

(2.34) 

(2.35) 


(2.36) 


Equation (2.32) implies that the potential may only jump by a 
constant across the shock. For convenience we shall assume that [cp] =0 
Note that u, v, w are the perturbation velocity components and are all 
zero upstream of the shock. On the shock surface the value of the 
perturbation potential is taken to be zero. 


2.3 FOEMULATION OF THE BOUNDARY VALUE PROBLEM IN CONICAL COORDINATES 


The flow field to be analyzed here is conical, that is, flow 
properties are constant along radial lines extending from the origin. 
For this type of a flow the equations derived in the previous section 
assume special forms. Since the potential is constant along rays it 
will be a function of the coordinates — and ^ . Thus we can define 

X X 


conical variables 







Y = -^ = -5- 

Bx bx 
where B = 

In terms of the above conical variables the perturbation potential 
assumes the form 

cp(xt y, z) = X 4)(Z, Y) (2.37) 

and the transonic equation (2.25) becomes 

ZZ ‘^YY^ 
(2.38) 


(2.39) 

(2.40) 

(2.41) 

It proves useful to write the equation for the perturbation 
potential in a curvilinear coordinate system. The elliptic cylinder 
coordinates, which were chosen, have the advantage of fitting well 
into the numerical scheme. 


> 11 ^ 


[K + (y+ 1) (4>-Z(t)2-Y(j)Y)[ JZ^<|)22 + 2 ZY(J>y2 + Y^(J)yy 


Since the wing is also conical 


f(x, -g) = xF(Z) 


and this makes the boundary condition on the body 


= B{F(Z) - ZF'(Z)} 


The two boundary conditions on the shock then are 

- K(Z<f)2+Y<f)Y)^- ^ (Z(|)2+Y<J)Y).^j+ ^ (c))^ + = 

' / B 


(j) = 0 
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The elliptic cylinder coordinates are given by 
Z = cosh ^ COST! 

Y = sinh ^ sintl 

for ^ s 0 and 0 s ri < 2 tt 

Since the wings to be considered are symmetric with respect to the y 

7T 

and z axes, the range of r} is 0 ^ r) ^ shown in Figure 2.3. The 

wing extends from zero to one on the Z axis = 0) and the shock is 
a curve in the plane; whose shape and location are part of the solution 
for the potential function. 

In the new coordinates the equations for the velocity potential 
and the boundary condition on the wing become 


k + (y+1) 


- y (sinh C cosh ^ <t>^ - sinr) cosq 

( 1/22 2 2 \ 

■j (sinh ? cosh 5 “ 2 sinh ? cosh ^ sinri cosq + sin q cos 


(7 


sinh C cosh C f (sinh^^ cosh^^ - sin^q cos^q) + 


1 - y (sinh^C + cosh^C) j <t>^ + 


') 


fr 


^ 2 2 2 2 
sinq cosq ( (sinh ^ cosh ^ - sin q cos q) + 


1 2 2 
1 + y (cos q - sin q) 






(2.42) 


where 

2 2 

J = sinh C + sin q 

and 

4>^(0,q) = ((j^Ccosq, 0) sinq (2.43) 
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On the shock surface the boundary conditions take the form 

1 2 yf 1 3 ) 

(sinh^ cosh^cj)^ - sinn cosr\<l)^) - (sinh^ cosh^4>^ - sinri cosncj)^) j 

+ — ^ ^(()^ + (f)^^^sinh^^ cos^n + cosh^^ sin^D^ = 0 (2.44) 

B K 

and (j) = 0 on the shock (2.45) 

Since the potential and the shock are sjnmnetric with respect to 
the Z and Y axis 


d) = 0 


at 


n = 0, ^ 


(2.46) 


0 

dn 


at 


n = 0 , 5 


(2.47) 


Fig. 2.4 gives a summary of the boundary value problem in the 
cross plane. 

Various physical quantities caii be derived from the solution of 
the transonic potential equation. Two which shall be considered here 
are the pressure coefficient Cp and the drag coefficient C^. 


The pressure coefficient is defined by 

p - Poo 2 

^ ~ 1 2 ~ 2 
4 p U ym 

2 ^00 ' oo 


fe-‘) 


(2.48) 


but 


2^ = 


Y 

2 \Y-1 


Poo \a^ 
\ 00 


and using equations (2.16) and (2.21) we obtain 


Cp 2 

YM 


K- 


(Y-l)Mf 6^/3 cpj - li 
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For small perturbations this expression may be expanded in series the 
result to lowest order being 

C = - 26^^^ (p 
P X 

In conical coordinates the pressure coefficient on the wing at x = 1 
takes the form 


Cp = - 26^/^ ((})- Z<^^) 


(2.49) 


The drag coefficient is obtained by integrating the pressure 
coefficient over the surface of the wing and normalizing to the area 
of the wing plan form. 




/A 


dA n • i 


(2.50) 


dA is an area element on the wing surface and n is the unit vector 
normal to the wing surface. In conical coordinates equation (2.50) 
becomes 




f 


Cp(F - ZF')dZ 


(2.51) 
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CHAPTER 3 


NUMERICAL PROCEDURES AND RESITT.TS 
3.1 INTRODUCTORY REMARKS 

The transonic equation for the perturbation potential (equation 
2.42) has the general form 

ai 4)^^ + 2a2 + 33 <t>^^ + = 0 (3.1) 

where the coefficients a^ through a^ are function of U» 4 >^» and 
the similarity parameters K and B. Depending on the sign of the 
quantity 

D = a^ - ai a 3 

equation (3.1) is elliptic parabolic or hyperbolic if D is negative, 
zero or positive respectively. 

For the cases to be considered here, equation (2.42) is of elliptic 
type within the region bounded by the shock wave (Fig. 2.4) and centered 
differences are proper to approximate its derivatives. 

The shock shape and location are not known a priori. They are 
part of the solution of equation (2.42). With this in mind a method 
has been developed that solves simultaneously for the shock shape and 
location, and the flow potential in an iterative manner. This is done 
by initially locating the shock near the Mach cone and then proceeding 
as follows: 

(i) Solve equation (2.42) subject to the boundary conditions on the 
Y and Z axes and <j) = 0 on the shock. 

(ii) Split equation (2.44) in such a v/ay that a new shock location is 
obtained using Information from step (i) 


(iii) Go back to step (i) with the shock at the new location and repeat 
the sequence, until the shock location converges. 

Equation (2.44) is split in a way that makes the sequence just presented 
a convergent one. 

The following sections give the details of the procedure just 
described. 


3.2 PRESENTATION OF DIFFERENCE FORMULAS AND RELAXATION ALGORITHM 


Consider the numerical solution of equation (2.42) with the shock 
near the ?fech cone. The location of the Mach cone is known from linear- 
ized theory; at x = 1 it is given by 

; 

. where 


On the Z axis 

. Z - cosh5„c 

and the initial shock location is taken to be 



cosh 


-1 


1 

Bv^ 


Two types of boundary conditions occur. On the shock the 
potential is specified and this is incorporated directly into the 
difference formulas by taking j ~ ^ other 

boundaries derivatives are specified and these are incorporated on mesh 
points which are offset one half mesh spacing from the actual boundary. 


This gives a simpler and more accurate treatment of the boundary 
condition. ^ 

The grid geometry is rectangular as shown in Figure 3.1 and the 
difference equations used for the interior points (i;^l, 3^1, jniax) 
are 






♦SEI.J 


4> 


nni.j 


" 

. .Vi.j - ^ , o(,e2^ 

l^K 

^i 1+1 ~ ^^i 1 *^i 1-1 2 

= ^ ±LL-L + 0(An ) 


(3.2) 

(3.3) 

(3.4) 

(3.5) 


'Cni.j 


+ O(AgAn) (3.6) 


where the subscripts l,j give the location where the solution is 
computed. These difference equations are second order accurate. 

For grid points along i = 1, j = 1, or j = jmax the difference 
formulas are modified to incorporate the specified values of and 
which are given along ^ = 0 and n = 0, -j > respectively. To Illustrate 
this consider a grid point at 1 = 1 for all j. This point is a distance 
away from the boundary ^ = 0 where is knovm from equation (2.43), 


and is denoted by PC(^). 
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With a fictitious mesh point at i = - 1 



c()^ at 1 = 1 can he written as 


*2,.1 - *2,.i -*• 

2A5 2AC 


* 2.1 - * 1 , 


+ 1 * 1,.1 - *- l . 

2 A? 


giving 1 


(3.7) 


is evaluated in a similar fashion, giving 


<J>o . - 


-PC(j) 0>2. 1 


TO.J 


PC(j) (3.8) 


For the mixed derivative at i = 1 three different forms appear 
depending on the value of j . If j 1 or j max we can make the approx 


imation 


and use equation (3.7) to obtain 


%i.j = is (pca«)-pca-i) + (3.5) 

If i = 1 and j = 1 (or jmax) , <()g^ is computed by averaging the values 
of at neighboring points as shown below. 


1 3 , 


,1 


I 


I 2 2 2 


'f ' 

I ^f’fl 


* — 


'll '2 1 

+ Wr i - 

I ^ ' , 


l.i 


2 . '2 




Cni.i 


= r/<^ 


+ <p 


+ <)) 


+ (f> 


m, T 1 ^Tii, 1 


but 


<1) -, = 0 since d) = 0 on h = 0 

cm. j " 


and 4)^^ on 5 = 0 can be calculated from equation (2.43) and will be 
denoted by PCA(j) giving 


<P 


cm.i 4 


= T PCA(l) + i PCA(i) 


. * 1.1 - * 2.1 - * 1 . 2 * * 2.2 

4A?An 


(3.10) 


Similarly special forms for and on j = 1 and j = jmax for' 

all i, and (|)^^ for i = 1, j = jmax are derived with the result summar- 
ized in Table 3.1. 

The difference equations just presented are used in equation (3.1) 
and a difference equation for <p^ ^ is obtained. This equation is 
linearized, by evaluating the coefficients a^, B 2 and a^ using values 
of J from the previous iteration, and solved by a point relaxation 


algorithm whose formula is 

+ (1 - co)4,<” 


(n) 

j 


,(n) 


where (fi; ^ is the value at the nth Iteration 

i.j 




(n+l)th 


^(n+1) relaxed value 

J 


0 ) is the relaxation parameter 

3.3 MOVING THE SHOCK 

On the shock 4> = 0 and we can write 


<p^ d^ + (|)^ dn = 0 


(3.11) 


or 




M = _ In 

dn <f)^ 


on the shock 

and upon substituting the above in equation (2.44) and rearranging it 


we obtain the relation 













slnh^C cosh^^ = -|2^^jsinhC coshC sinri cosri + sin^ri cos^ri 


(t>i- fsinh? cosh^ + sinri cosri 


2JK '•'C 


WJ 


B^K 


^sinh^C cos^ri + cosh^^ sln^rijj 


(3.12) 


With the solution just obtained the right hand side of equation 
(3.12) is computed and a corrected new shock location is obtained 
whose shape is usually non elliptic. Since the same number of mesh 
points are kept, so that the shock is always on grid points, the new 
grid geometry is non-rectangular and is now a function of r|. A 
subscript J will be used, to show this dependence. Note that grid 

points which have the same value for the subscript i are not aligned on 
vertical lines any more. 

The values for the potential (with the shock at ^^) are 
transferred to the new grid points and are used as starting values fpr 
the relaxation procedure. 

In calculating r| derivatives (at a given i) values of (j) on a 
vertical line are needed, thus requiring interpolated values of (J> at 
points where the solution is not computed. Figure 3.2 shows these 
points (marked by x's) when the derivative is calculated at i,j . 

The values of the potential at these points are calculated applying 
a second order accurate interpolation formula which uses values of the 
potential at the three nearest grid points. For example the value of 
(() at a point s away from i, j+1 is 


29 



♦s,J+l ” ® (♦l.J+l ■ ^ 2 ' ('•'l+l.J+l ' 

This interpolation formula has the same accuracy as the difference 
equations. The solution for (|) Is analytical in the neighborhood of the 
point i,j+l, therefore the overall accuracy of the numerical scheme is 
preserved. 

3.4 COMPUTATIONAL DETAILS 

Sections 3.2 and 3.3 give the basic method for solving the 
transonic potential equation. In applying these principles certain 
points must be resolved. These include choice of relaxation parameter, 
number of grid points used as well as the determination of convergence 
of the iterative scheme and overall accuracy. 

The choice of the relaxation parameter (co) determines the rate of 
convergence of the algorithm. For Laplace’s equation an optimal value 
can be found from the spectral radius of the Gauss-Seidell method. 

The transonic equation is nonlinear and such a value must be found by 
experiment since it does not have an explicit form. The value of u 
that gave the fastest convergence (fewest number of iterations) for a 
grid with 40x20 points, in the ? and u respectively, B = .2 and K = 3.0 
is 1.93. This value was used, for all the computation since co did not 
depend strongly on B and K. 

The convergence of the solution, for a particular shock location, 
is based on the value of the residue - the difference in the potential 
at a given grid point on two successive iterations. . When the value of 
the residue at all points in the field is less than 2.5 x 10 ^ the 
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solution is assumed to have converged. The same criterion is used 
for the shock location with 10 ^ as the reference tolerance. Table 3.2 
gives the number of iterations needed for convergence (n) and the shock 
location (value of -^) on n = 0. The location prescribed initially to 
the shock wave is 2.65000; the rest are corrected locations obtained 
from equation (3.12). 


Table 3.2 

Successive Shock Wave Locations 


K = 3.00 

B = .2 


z 

n 

b 

118 

2.65000 

80 

2.79355 

66 

2.86181 

62 

2.89161 

48 

2.90411 

21 

2.90925 

9 

2.91132 

5. 

2.91216 

2 

2.91249 

1 

2.91265 


2.91267 


The accuracy of the methods has been tested by mesh variation. 

For shock locations more than one span away from the tip of the wing 
the grid used has 40 points in the ^ direction and 20 in the r) direc- 
tion. However, for cases with the shock closer to the tip of the wing 
a 20 X 20 grid is sufficient. 



# 




f 


# 










In calculating the pressure coefficient, values of the potential 
on the wing are needed. These are evaluated by extending the solution 
to the wing surface using a second order accurate extrapolation equation 
which has the formula 



The drag coefficient integral, equation (2.51), is evaluated numerically 
using the trapezoidal rule. 

This section shall conclude with a word about computation times. 

Two different time scales appear in this problem. Consider starting a 
computation for one set of conditions (B and K) from the solution to a 
different problem. If in the new problem the shock location is signifi- 
cantly different a well converged solution may take 15 cycles (shock 
location corrections.) If solutions are sequenced so that shock loca- 
tions do not change by much the same degree of convergence may take 
only 8-10 cycles. A computation that requires 10 cycles, on a grid of 
40 X 20, takes about 40 seconds on an IBM 360/91 using double precision 
arithmetic. 

3.5, RESULTS 

In this section computed results based on the theory of the 
previous sections is presented. Two different conical wings, one of 
rhombic cross section the other circular, are considered for various 
values of the similarity parameters K and B. The methods can be used 
to compute other wings which are pointed at the tip . 

In conical variables the rhombic wing is given by 

F(Z) = 1- Z 0<Z<1 (3.13) 
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(3.14) 


This equation corresponds to 
y = 6xF(Z) = fix 

in the physical coordinates and simplifies to 
y = « (l - 
at X = 1. 

For the circular wing 

F(Z) =1-Z^ 0<Z<1 (3.15) 

and in the physical coordinates it becomes 



at X = 1. 

Figures 3.3 and 3.4 give the shape and location of the shock wave in 
in the cross plane. Figures 3.5, 3.6 and 3.7 give the pressure coeffi- 
cient along the wing surface for the indicated values of K and B. 




Figure 3.3. Shock Waves for a Rhombic Wing at K = 2.0. 
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K = 3.0 
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of a Rhombic Wing for K = 1.0. 
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Figure 3.6. Pressure Distributions on the Surface of a Rhombic Wing for K = 2.0. 
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Figure 3.7. Pressure Distributions on the Surface of a Circular Wing for K = 3.0. 
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CHAPTER 4 


SLENDER WINGS 

4.1 SLENDER WING THEORY 

In this section the theory of "slender" wings will be presented. 

This theory has been studied by several authors, but the formulation 
12 

given by Messiter for the non-lifting case is particularly suitable 
for our purpose. The problem can be defined in several ways, all of 
which are equivalent. From a mathematical point of view the fundamental 
assumption is that the solution for the potential near the wing (inner 
expansion) is essentially different from the solution near the shock 
wave (outer expansion) . Using a physical picture one could take the 
reduced aspect ratio, b/n^ - 1, as a basic parameter and discuss the 
solution as it tends to zero. Another possibility is to start from non- 
slender wings and require the parameter B to approach zero with 6, 
subject to the restriction that ^ also approaches zero in order that 
the wing remain thin. The last choice is the most convenient and is 
used to obtain inner and outer expansions for the potential. The inner 
expansion corresponds to the solution of the Laplace equation in the 
cross plane. Incompressible flow, and the outer expansion to axisym- 
metric, non-linear transonic, flow about an equivalent body of revolu- 
tion. 

This equivalent body is defined as the body of circular cross- 
section which has the same longitudinal distribution of cross-sectional 
area as the wing. The equivalent body is defined by 

r =6^ E(x) (4.1) 

PREX)BDING page blank not filmed 
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where 


,/~ 2 ~ 2 

r = Vy + z 

E(x) is the longitudinal cross-sectional area distribution and 6 is 

e 

the equivalent radius at x = 1. 

The cross-sectional area of the wing, at x = 1, expressed in terms 
of the equivalent radius is related to the wing dimensions by the con- 


stant of proportionality k 
•n-6^ = 7Tkb6 = Trk6^'^^ B 


(4.2) 


For the outer expansion 

(p (x,r ; K^) = 6^ |e(x)E'(x) logr* + g^^(x; (4.3) 

r* = r6 

e 

- 1 

K = — (similarity parameter for axisymmetric flow) (4.4) 

6 r ^ 


and for the inner expansion 


(^(x,y,z; K ) = $1 cp log 

e el e e z 


(4.5) 


where 


— y — z 

y b * " = b 


cp^= E(x)E’(x) 


(4.6) 




_ fZ =— Z- 

= ^ / k^(x,cr)logVy^ + (z - a)^ da + g^(x; Ke) 


(4.7) 


*s^(x) 


gj^(x; K^) is the essential unknoTO part of the pressure distribution 
and can only be determined from the solution of the outer problem. The 


# 


# 


2 

limits of the integral describe the wing plan form and y = 6h(x,r-) 

D 

gives the shape of the wing. 

The pressure coefficient Cp is obtained from the inner solution 
12 

and has the formula 

P - - 2 106(b«^) + 

e 

4.2 CONICAL SLENDER WINGS 


The equations for slender wings presented in the previous section 
can be simplified using the properties of conical flows. For conical 
wings . 

E(x) = X 

h(x,a) = xFO:) , ^ ^ 

h^(x,a) = F(E) - E F’(E) 

and since the potential is a conical variable 

gj^(x; K^) = - xlogx + xCj^(K^) (4.9) 

With the above information, equation (4.7) can be written as 


^2 


X 

iTk 



(F - EF') 


I log X 


+ log 



+ (z 



- X log X + xC^(K^) 


(4.10) 
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and since 


j. 

i- f 

TTk / 
-1 


(F - ZF’)dZ = 


TTk 


1 

/ 


FdE = 1 


1 

^2 = ^ f (F - log |Z-^I 


(4.11) 


on Y = 0. 


2 2 /'^ 

From equation (4.2) 6^ = kB(S so 

K = ^ 

^e kB 


(4.12) 


and 


b6 = 
e 


The expression for the pressure coefficient at x = 1 then becomes 

s 

62 " ■ 
e 


jlog k + 3 log B + 2 V? 2 x| 


or in terms of 6 


= - kB log k + 3 log B + 2<P 


I 


2x 


(4.13) 


For the rhombus k = — and F = 1 - Z thus giving 

^/3 " ■ f ® f + 3 log B + log(l-z2) - 2 + 2C^(K^)| (4.14) 

8 2 

The result for the circular arc where k = and F = 1 - Z is 


1/3 = - If ® If ^ 3 log B + z2 - 1 + 


I log (l-z2) - I + 2C^(Kg,j 


)! 


(4.15) 
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The drag coefficient can now be calculated using the above two 
equations and performing the integration in equation (2.51) gives: 
Rhombus 

^/3 " f ® 2 log 2+3 log b| - 2C^(K^)| (4.16) 

Circular Arc 

^/3 - if ®{ir - T If + 3 log 2 + 3 log b) - I C^(K^)j(4.17) 

In order to make the formulation for the pressure and drag 

coefficients complete, we must compute C. (K ). This can be done in 

i e 

several ways. One method of obtaining Cj^(K^) is to Integrate the first 

14 

order transonic equation for flow past a slender cone. A second 

method is to use the value of the pressure coefficient, on the surface 

13 

of a slender cone from Kopal's table, in the similarity law given by 
the equation^^ 

4 

- 2C (K ) = -r + 4 log a - 1 

® a 

Cp is the pressure coefficient on the surface of the cone, and a is the 
tangent of the cone semi-vertex angle. 

Yet another method is to compare the pressure coefficient obtained 
ntnnerically, for a particular shape, with the slender wing formula for 
that particular shape and choose C]^(Kg) such that the two have the same 
value at a given location on the wing. The location that was chosen is 
at Z = 0. This method also serves as a check on slender body theory 
since the results obtained numerically represent the solution of a 
higher order equation. 


In Figure 4.1, 2C^(K^) computed using Kopal's table is indicated 
by the solid curve whereas the results obtained by matching are shown 
by triangles and squares . 

4.3 COMPARISON OF RESULTS 

The comparison between the results obtained numerically and the 
slender wing solution, equations (4.14), (4.15), (4.16) and (4.17), is 
shown in Figures 4.2 - 4.7. 
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Figure 4.1. The Function C, (K ). 
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Figure 4.3. Pressure Distributions on the Surface of a Rhombic Wing for K = 2.0. 
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Figure 4.6. Pressure Distributions on the Surface of a Circular Wing for B = .2. 
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Figure 4.7. Drag Coefficients for Rhombic and Circular Wings for B = .2. 
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CHAPTER 5 


CONCLTJBING REMARKS 

The method presented in this dissertation falls in the general 
category of shock fitting techniques. It is simpler than conventional 
methods mainly because the state on one side of the shock is known. 

The results of the computations indicate that, generally the method has 
a wide range of applicability. The limitation of the algorithm is best 
expressed in terms of the location of the shock wave, i.e. it must not 
fall on the wing. 

For wings which have rounded tips, such as conical wings of 
elliptic cross section, the same algorithm can be used provided that the 
numerical solution is modified in the neighborhood of the tip in order 
to recover the correct singularity. The conical solution in the neigh-^ 
borhood of a rounded tip is given by the solution, in the hodograph 
plane, of two dimensional transonic flow towards that tip. 

The agreement between slender body theory and the numerical cal- 
culations indicates the usefulness of equations (4.14) - (4.17), (II-5) 
and (II-6) . These equations are applicable for values of for which 
Cj^(K^) has a similarity form. Expressed in terms of the similarity 
parameters B and K, C^(K^) possesses such a form if B < .20 and K < 2.0. 

It should be noted that the wing shapes considered are given by 
the family of curves F(Z) = 1 - Z^, n = 1, 2, 3, 4 which have dif- 
ferent area distributions along the Z axis. Also, the function 
C^(K^) does not seem to be particularly partial to any of these 

shapes . 
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The computer program has been written to ease the understanding 
of the logic rather than optimizing it for computing speed. 

There are many areas in which future work would be beneficial. 
One V 70 uld he to extend this method to lifting wings. In particular 
conical xi/ings with arbitrary afterbodies could be studied. This can 
done by extending the solution from the conical region to the region 
over the afterbody using the method of characteristics. 
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APPENDIX I 


PROGRAM LISTING 

The computer program used in this dissertation was written in 
Fortran IV for an IBM 360/91. The main program is included here for 
reference. 

The following subprograms are called by the main program: 

INTERP Computes interpolated values for the potential 
OUTPUT Prints out values of the potential 

SLEND Computes values of the pressure coefficient using slender wing 


theory. 
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IF( I .EQ. 1) GO TO ?00 
I F( J, EO. 1 ) GO TO 230 
IF< J,EQ. JMAX) GO TO 240 
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702 DO 703 J=1.JMAX 

IF(J,EQ«1) C A( J )=,5*ODA* ( CMAX( J +l #L )-CMAX(J,L 
IF(J.EQ.l) GO TO 705 



1000 CONTINUE 
3000 CONTINUE 

WRITE(6,709)L 





WRITE ( 6. 1 00 



900 CONTINUE 

WRITEI6. 1002) 

WRITEI6.899I SUM 

899 F0RMAT(1H0,' DRAG = *,F15.7) 








CALL SLEND(TK,B,CPM) 


601 WRITE(6,602) 

602 FORMAT(» THE SOLUTION DIVERGED*) 

603 CONTINUE 

OOP CONTTNUF 

STOP 

END 
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APPENDIX II 

In this appendix two more nonlifting conical wings will be considered 
and the function will be calculated using the matching method 

described in Section 4.2. 

In terms of conical variables the first is described by the cubic func- 
tion 

F(Z) =1-Z^ 0<Z<1 (II. 1) 

which in the physical coordinates corresponds to 
/ 3 \ 


y = 6x ( 1 - 


3 

b X 


(II. 2) 


and the second by the quartic function 
F(Z) = 1 - Z^ 
which corresponds to 


0 < Z < 1 


(II. 3) 


y = 6x { 1 - 


b X > 


(II. 4) 


Formulas for pressure distributions on slender wings can now be 
obtained using equations (4.2), (4.8) and (4.11). The constant of pro- 

3 

portionality k, obtained using equation (4.2), is equal to — for the cubic 
and for the quartic. Performing the integration in equation (4.11) and 
the differentiation in equation (4.8) gives 




log + 3 log B + -| Z + 2Z^ + 4 log (1-Z^) 




3 

( I~2) _ A + 2C TK ) 
(1+Z) 3 + 2Cj_(K^) 


for the cubic and 
C 


p _ 16 

j2/3 5 tt 


B 


log 41 + 3 log B + 41 + 3Z^ + log (1-Z^) 

(1-Z) 


+ 1 tSt - i 


for the quartic. 


(II. 5) 


(II. 6) 
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The values of the function Cj^(K^) obtained by matching equations (II. 5) 
and (II. 6) with the pressure coefficients obtained numerically are shown in 


Figure 4.1. 


